General Relativistic Decompression of Binary Neutron Stars During Dynamic Inspiral 
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We investigate the dynamic stability of inspiraling neutron stars by performing multiple-orbit 
numerical relativity simulations of the binary neutron star inspiral process. By introducing eccen- 
tricities in the orbits of the neutron stars, significant changes in orbital separation are obtained 
within orbital timescales. We find that as the binary system evolves from apastron to periastron 
(as the binary separation decreases), the central rest mass density of each star decreases, thus sta- 
bilizing the stars against individual prompt collapse. As the binary system evolves from periastron 
to apastron, the central rest mass density increases; the neutron stars re-compress as the binary 
separation increases. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.40.Dg, 02.60.Cb 



INTRODUCTION 

Detailed and accurate models of binary neutron star 
coalescence phenomena, from quasi-equilibrium orbits 
through plunge and merger to the subsequent formation 
of the final compact object, will be required in order to 
extract information regarding the structure of neutron 
stars from detected gravitational wave signals. One de- 
tail of some debate over the past decade has been the so- 
called "neutron star crushing effect" , reported in P, 0, [1] . 
This effect was first reported in where binary neu- 
tron star simulations employing a variety of simplifying 
assumptions indicated that as the binary stars spiraled 
inwards, the general relativistic gravitational interaction 
between the two stars caused a destabilization, trigger- 
ing the collapse of each star to individual black holes well 
before the plunge and merger phase of the inspiral. Since 
then, studies using several different sets of approxima- 
tions 0, H 0, 0, S, S 0, EH predict the exact oppo- 
site, namely that the gravitational interaction between 
the two stars act to stabilize each star as the binary sep- 
aration decreases during inspiral. An error in the formu- 
lation used in the original "neutron star crushing effect" 
studies was pointed out in [l^ . However, subsequent 
studies still claim a destabilization effect 

In this Article, we report on the first fully dynamical 
general relativistic simulation results aimed at studying 
this "neutron star crushing effect" for inspiraling binary 
neutron stars. In the past several years, fully general 
multiple-orbit binary neutron star simulations have been 
performed in numerical relativity [l^ [13. We use the 
evolution code previously described in [13| to perform 
multiple-orbit simulations of binary neutron stars (for a 
detailed analysis of the accuracy of the simulations per- 
formed in this Article, see [IH|). We begin each simula- 
tion with initial data corresponding to quasi-equilibrium, 
circular orbit configurations modified such that the initial 
orbital angular velocity of the neutron stars is decreased 
by varying amounts. Each initial data configuration sat- 
isfies the initial value constraints of general relativity. 



The simulations resulting from these eccentric-orbit ini- 
tial data sets permit the study, within the context of full 
general relativity, of the compression/decompression ef- 
fect by correlating the central rest mass density of the 
neutron stars with the proper separation of the stars as 
the evolution progresses through several apastron (max- 
imum separation) / periastron (minimum separation) 
points during the inspiral. We find that the stars do, 
in fact, stabilize as the binary separation decreases. This 
result is the first fully dynamic general relativistic demon- 
stration of the decompression of binary neutron stars; 
the typical simplifying assumptions made in the treat- 
ment of the problem to date, such as the truncation 
of a post-Newtonian expansion, the imposition of quasi- 
equilibrium conditions, or the forcing of the spin states 
of the evolution of the neutron stars (e.g., forcing coro- 
tation or irrotation during the inspiral) , are not made in 
the analysis presented here. 



SIMULATION RESULTS 

The numerical relativity / general relativistic hydrody- 
namics code used here is described in [l3| and analyzed 
for accuracy in ^15| . The grid size used for simulations is 
643 X 643 X 325, employing a mirror symmetry about the 
equatorial plane of the neutron stars, z — Q. We choose 
units such that the gravitational constant G and the 
speed of light c are identically 1. A spatial discretization 
of Ax/m = 0.148 is used, which corresponds to roughly 
40 points across the diameter of each neutron star. The 
boundary of the domain is thus eight neutron star diam- 
eters away from the center of mass of the system. The 
physical parameters of three simulations NS-1, NS-2, and 
NS-3 are shown in Table HI The initial data corresponds 
to the quasi-equilibrium, corotating, circular orbit initial 
data described in except that here, the initial orbital 
angular velocity parameter Oq has been decreased from 
the circular orbit value by a factor of 1%, 2%, and 3% for 
simulation NS-1, NS-2, and NS-3, respectively. The re- 



Config. 




T / 2 


Mo/m 




NS-1 


0.01618 


1.121 


1.073 


18.10 


NS-2 


0.01599 


1.094 


1.073 


18.12 


NS-3 


0.01580 


1.067 


1.073 


18.14 



TABLE I: The physical parameters of the initial data used for 
numerical relativity simulations NS-1, NS-2, and NS-3. These 
initial data sets correspond to quasi-equilibrium, corotating, 
circular orbit initial data with the initial angular velocity fio 
reduced from the circular orbit value by 1%, 2%, and 3%, 
respectively. The symbols Jo and Mq are used to denote the 
initial angular momentum and total rest mass of the con- 
figuration, respectively. The initial proper separation of the 
neutron stars, defined by Eq. 57 in [131, is denoted as (rp)g. 
A polytropic equation of state with adiabatic index F = 2 
is used. The rest mass of each neutron star corresponds to 
82% of the maximum stable TOV rest mass configuration. 
All physical quantities in this Article are expressed in units 
of m, which we define to be twice the ADM mass of a single 
stationary neutron star configuration with rest mass Mo/2 
uniformly rotating with angular velocity Qo- 

suiting initial data satisfies the Einstein equation initial 
value constraints, and results in binary evolutions whose 
orbits contain increasing amounts of eccentricity. 

The initial spin state of each of the neutron stars cor- 
responds to a corotation with the initial orbital angular 
velocity. For binary neutron stars of 1.4 solar masses 
each, this corresponds to a spin period of approximately 
5 milliseconds. However, it is important to note that 
during the dynamical evolution, the spins of the neutron 
stars do not remain tidally locked to the orbital angular 
velocity (see Fig.[T]). Previous quasi-equilibrium sequence 
inspiral studies have assumed either strictly corotating or 
strictly irrotational neutron star spin scenarios; in these 
cases, the spin states of the neutron stars are fixed to 
have either corotating spin or to be irrotational during 
the inspiral sequence. By numerically solving the general 
coupled Einstein/hydrodynamics equations, we are able 
to study the more astrophysically relevant case where 
the non-zero spin of each neutron star actually evolves in 
time via the equations of motion. 

In Fig. [21 the angular momentum of the binary sys- 
tem simulations is plotted as a function of coordinate 
time; tcoord = 500 m corresponds to roughly two orbital 
periods. For comparison, the angular momentum for so- 
lutions to the 2.5 post-Newtonian (dotted curves) and 
4.5 post-Newtonian (dashed curves) point particle equa- 
tions of motion is plotted along side the simulation re- 
sults (solid curves). The overall rate of decrease in angu- 
lar momentum for the numerical relativity simulations is 
roughly equivalent, if not slightly more, than predicted 
by the post-Newtonian equations of motion. However, 
the difference in total angular momentum between the 
numerical relativity simulations NS-1, NS-2, and NS-3 
(solid curves) is nearly twice as large as the difference be- 
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FIG. 1: The circulation C = §^ da hu^,{-^f (h is the rel- 

ativistic specific enthalpy of the neutron star fluid, u'^ is the 
4- velocity of the neutron star fluid, and ct is a Lagrange pa- 
rameter labeling points on the curve) along constant rest mass 
density curves Cp in the equatorial plane of the neutron stars 
is plotted as a function of rest mass density for three dif- 
ferent neutron star binary configurations. The first (solid) 
curve shows the circulation of the initial corotating, quasi- 
equilibrium binary configuration of simulation NS-3, while 
the second (dotted) curve shows the circulation of the binary 
configuration of simulation NS-3 at the time corresponding to 
the first periastron point. For comparison, the circulation of 
a corotating, quasi-equilibrium, equivalent rest mass binary 
configuration with a proper separation equal to that of simu- 
lation NS-3 at periastron is plotted (dashed curve). Clearly, 
the spin of the neutron stars during the numerical evolution 
NS-3 remains roughly constant. This is in stark contrast to 
a corotating quasi-equilibrium sequence in which the neutron 
stars are artificially "spun-up" as the binary separation de- 
creases. 

tween the post-Newtonian curves, even though the post- 
Newtonian initial data is selected using the same method 
as in the numerical relativity case, namely, using circular 
orbit initial data with a decrease in initial angular veloc- 
ity of 1%, 2%, and 3% (see Table This is most hkely 
due to the nonlinear effects of resolving the constraints of 
general relativity after decreasing the initial angular ve- 
locity parameter Hq when preparing the initial data for 
the numerical relativity simulations. 

Contour plots of the rest mass density for the initial 
data, as well as the first periastron points of simulations 
NS-1, NS-2, and NS-3, are shown in Fig. [3 The first pe- 
riastron points occur just slightly before the completion 
of three-quarters of an orbit for all three simulations. 



DECOMPRESSION OF BINARY NEUTRON 
STARS 

In Fig. [31 we plot the proper binary separation Vp (de- 
fined to be the spatial geodesic distance between the max- 
imum rest mass density point of each neutron star, see 
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FIG. 2: The total angular momentum J is plotted for the 
binary neutron star numerical relativity simulations NS-1, 
NS-2, and NS-3 (solid curves). All simulations correspond 
to roughly two orbital periods. For comparison, solutions to 
the post-Newtonian point particle equations of motion are 
shown, using the same initial coordinate separation, mass, 
and initial orbital angular velocity (i.e., 1%, 2%, and 3% 
below the circular orbit angular velocity) as the numerical 
relativity simulations. Dotted curves are solutions to the 
2.5 post-Newtonian (quadrupole level) equations of motion 
(see, e.g., |_L0|). Dashed curves are solutions to the 4.5 post- 
Newtonian equations of motion, which include 3.5 pN [l^ 
and 4.5 pN [13] radiation reaction terms, along with the 3.0 
pN conservative terms. 
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FIG. 3: Contour plots of the rest mass density p in the x — y 
equatorial plane of the binary neutron stars is shown for initial 
data (upper left) and the first periastron point for simulations 
NS-1 (upper right), NS-2 (lower left), and NS-3 (lower right). 
The central contour corresponds to a value of p = 0.0083/m^. 
Six more contours are shown; each successive contour corre- 
sponds to a decrease in rest mass density by a factor of two. 



Eq. 57 in |13i] ) as a function of the proper time tp mea- 
sured by observers located at the maximum rest mass 
density points of the neutron stars. An evolved proper 
time of tp — 300 m corresponds to roughly two orbital 
periods. Also in Fig. [4] is a plot of the maximum rest 
mass density Pmax of the neutron stars as a function of 
proper time tp. Oscillations in Pmax{tp) corresponding 
to the fundamental radial oscillation mode of each neu- 
tron star have been filtered out (these oscillations have 
a proper period of Tp = 35.0 m). Fig. 2] displays a clear 
correlation between the maximum rest mass density of 
the stars and the proper separation of the stars; smaller 
proper separation Tp corresponds to smaller maximum 
rest mass density pmax ■ More quantitatively, the correla- 
tion Cor[rp^ Pmax) between the proper binary separation 
and the maximum rest mass density is computed to be 

Qor{rp,pmax) > 0.95 (1) 

for all simulations. Convergence tests on simula- 
tion NS-3 yields a Richardson extrapolation value of 
Cor{rp,pjnax) = 1-00 (see Fig.E]). 

Fig. [3] also shows that the correlation between orbital 
separation and maximum rest mass density exists over 
gravitational radiation reaction timescales. This is evi- 
denced by comparing both the proper separation of the 
binary and the maximum rest mass density at successive 
apastron (local maximum in separation) points in Fig. U) 
First, we see that the proper separation Vp at the sec- 



ond apastron point (which occurs at roughly tp — 250 m 
for all three simulations) is approximately 10% smaller 
than the proper separation at the first apastron point 
(which occurs at the initial time tp = 0). Gravitational 
radiation emission drains the binding energy of the bi- 
nary, causing an overall decrease in separation from one 
apastron point to the next. We also see from Fig. 3] 
that the maximum rest mass density of the neutron stars 
is smaller at the second apastron point as compared to 
the first apastron point. Thus, these simulations also 
demonstrate the correlation between orbital separation 
and maximum rest mass density on gravitational radia- 
tion reaction timescales. 

In a post-Newtonian matched asymptotic expan- 
sion technique is used to obtain an expression for the 
fractional change in the central density pc of inspiraling 
binary neutron stars, perturbative in powers of the "tidal 
expansion parameter" a = R/r where R is the stellar ra- 
dius and r is the binary separation, such that 



as a — > 0. That is, when the fractional change in the 
central rest mass density of the stars is expanded about 
a = (i.e. infinite binary separation) in powers of a, 
post-Newtonian theory predicts the first non-zero term 
to be a^. We note that the limited range of binary 
separations, along with the large values of a obtained 
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FIG. 4: Top panel: the proper separation of the binary rp 
is plotted as a function of proper time tp as measured by ob- 
servers at the maximum rest mass density pmax of the neutron 
stars, for numerical relativity simulations NS-1 (solid curve), 
NS-2 (dotted curve), and NS-3 (dashed curve). Bottom panel: 
the maximum rest mass density pmax of the neutron stars is 
plotted as a function of tp for numerical relativity simulations 
NS-1 (solid curve), NS-2 (dotted curve), and NS-3 (dashed 
curve). Oscillations in pmax{tp) corresponding to the funda- 
mental radial oscillation mode of each neutron star have been 
filtered out. 

in our simulations {a = Rp/vp « 0.4, where Rp is the 
proper stellar radius computed from the proper stellar 

1/3 

volume Vp according to Rp = {ZVp/{'i'K)) ' ), prevent us 
from accurately determining the power of the first non- 
zero term in an expansion of the fractional change in 
the central rest mass density as in, e.g., Eq. [21 Simula- 
tions involving much smaller values of a (larger binary 
separations, which imply longer orbital periods), which 
will only be capable with adaptive mesh refinement tech- 
nology and/or higher order methods, will be required to 
verify the post-Newtonian prediction Eq. [2l 
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